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I.  INTRODUCTION 


When  designing  an  artillery  shell,  it  is  necessary  to  develop  a 
vehicle  which  will  fly  with  stability  under  a wide  variety  of  aero- 
dynamic conditions.  A range  of  propellant  charges  may  be  used  giving 
the  shell  launch  velocities  covering  a spectrum  from  subsonic  to  super- 
sonic. The  shell  will  also  slow  in  flight,  particularly  near  the  apex 
of  its  trajectory.  It  is,  therefore,  important  that  the  shell  fly  with 
stability  in  subsonic,  transonic,  and  supersonic  flight  regimes. 

Difficulties  are  often  experienced  by  projectiles  at  transonic 
velocities.  Aerodynamic  properties  such  as  drag  and  the  pitching 
moment,  which  is  critical  to  stability,  will  reach  peak  values  at  some 
transonic  Mach  number.  This  peak  can  form  in  a Mach  number  range  which 
may  be  limited,  for  example,  between  92  < M < .94.  The  sharpness  of 
this  critical  behavior  as  well  as  the  value  of  the  critical  Mach  number 
are  very  sensitive  to  body  geometry.  A slight  change  in  boattail 
length  may  make  the  difference  between  a successful  shell  and  one  whose 
behavior  is  unpredictable. 

Aerodynamic  range  and  wind  tunnel  testing  are  difficult  and 
expensive,  particularly  at  transonic  velocities.  Therefore,  it  is  of 
great  importance  for  artillery  projectile  design  to  develop  a computa- 
tional capability  which  can  provide  guidance  in  the  design  of  shell 
configurations  and  reduce  aerodynamic  testing  requirements.  Techniques 
have  been  established1  and  computers  are  now  available  which  should 
make  possible  the  development  of  useful  computational  design  tools, 
particularly  for  the  limited  and  seemingly  simple  geometries  found  in 
artillery  projectile  shapes.  The  techniques  presented  here  will  treat 
inviscid  transonic  flow  over  projectiles  allowing  the  computation  of 
lift  and  pitching  moment.  Magnus  effects  will  be  included  in  following 
work  when  viscous  effects  are  included. 

The  simplicity  of  the  shape  of  an  artillery  shell  is  somewhat 
deceptive.  There  are  sharp  corners  and  other  discontinuities  in  the 
surface  slope  and  curvature.  These  discontinuities  create  strong  shock 
patterns  at  transonic  velocities  as  may  be  seen  in  the  shadowgraph 
presented  in  Figure  1.  An  understanding  of  these  shock  patterns  is 
critical  to  an  understanding  of  the  aerodynamics  of  the  shell.  The 
shock  on  the  windward  surface  of  the  boattail  is  a great  deal  farther 
aft  than  the  shock  on  the  lee  surface.  This  pattern  generates  a strong 
upward  force  on  the  boattail  tending  to  overturn  the  shell  and  thus 
generates  a large  pitching  moment. 


1.  F.  R.  Bailey  and  W.  F.  Ballhaus,  "Comparisons  of  Computed  and 

Experimental  Pressures  for  Transonic  Flow  About  Isolated  Wings  and 
Wing  Fuselage  Configurations,"  NASA  SP-347 , Vol.  2,  March  1975, 
pp.  1213-1226. 
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The  resolution  of  the  shock  patterns  associated  with  flow  over 
the  corners,  makes  necessary  the  application  of  very  fine  grids  in  the 
numerical  procedures  used  to  obtain  predictions  of  the  flow  over  the 
shell.  The  necessity  of  making  three-dimensional  computations  with 
these  fine  grids  leads  to  the  problem  of  a large  computer  storage 
requirement  and  to  the  problem  of  a long  run  time.  Grids  on  the  order 
of  ten  thousand  total  points  have  been  needed  in  two  dimensions  to  suit- 
ably resolve  the  shock  patterns  and  corner  flows  associated  with  shell. 
Computer  time  for  such  a solution  based  on  a line  relaxation  scheme  is 
close  to  one  minute  on  a CDC  7600.  A carry-over  of  the  two-dimensional 
solution  to  three  dimensions  using  line  relaxation  techniques  could 
multiply  time  and  storage  requirements  by  a factor  of  30.  Such  a compu- 
tation time  would  be  unacceptably  long  for  engineering  design  purposes. 
Techniques  suited  to  cylindrical  problems  in  three  dimensions  will  be 
developed  in  this  paper  which  require  about  four  times  the  effort 
necessary  for  a two-dimensional  solution.  This  effort  leads  to  a run 
time  of  about  three  minutes  on  the  7600  which  is  considered  to  be  suit- 
able for  engineering  use. 


II.  THE  TRANSONIC  SMALL  DISTURBANCE  EQUATION 

The  approach  used  here  to  compute  transonic  flow  is  based  on  the 
solution  of  the  transonic  small  disturbance  equation, 

(1  - M2  - M2 (y*1)$2]$zz  ♦ $rr  ♦ 4>r/r  + 4>00/r2  * 0 (1) 


which  may  be  derived  from  the  Euler  equations  applied  to  a slender  body 
with  a pointed  nose2.  The  solution  of  this  equation  is  the  potential  4> 
from  which  velocities  and  pressures  in  the  flow  field  may  be  obtained. 
This  is  a non-linear  partial  differential  equation  of  mixed  elliptic- 
hyperbolic  type  written  in  a cylindrical  coordinate  system  (z,r,9)  as 
shown  in  Figure  2.  The  free  stream  Mach  number  is  given  by  M in  this 
equation  and  the  ratio  of  specific  heats  (1.4  for  air)  is  represented 
by  y . Boundary  conditions  are  given  on  the  surface  and  at  infinity  by: 


and 


surface 


= dR/dz 


4^  = ctrcos(6) 


In  these  equations  R is  the  radius  of  the  surface  and  a is  the  angle  of 
attack  of  the  projectile. 


2.  Holt  Ashley  and  Marten  Landahl,  Aerodynamics  of  Wings  and  Bodies , 
Addison-Wesley  Publishing,  Inc.,  Reading,  MA  (1965). 
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In  practice  it  is  more  convenient  to  use  a set  of  boundary 
conditions  in  which  all  the  important  information  appears  in  the  body 
surface  equation,  thus  leaving  the  boundary  condition  at  infinity  as, 

♦.  ■ 0 • 

This  desired  arrangement  of  the  boundary  conditions  may  be  accomplished 
by  the  transformation, 

♦ ■ v ♦ a r cos  (0)  , 
v • ♦ - a r cos  (0) 

This  transformation  does  not  alter  the  differential  equation  that  must 
be  solved.  It  does  alter  the  boundary  conditions,  however,  producing 
the  new  set: 


surface 


dR/dz  - a cos  (0) 


The  solution  to  equation  (1)  has  been  shown  by  Bailey3  to  give 
good  results  for  the  slender  body  case  at  zero  angle  of  attack.  This 
equation  has  also  been  studied  both  numerically  and  analytically  for 
many  years  and  it  is  simple  enough  that  much  valuable  insight  may  be 
gained  from  it. 


III.  FOURIER  TRANSFORMATIONS 

An  examination  of  the  closely  related  linear  equation  produced  by 
holding  the  coefficient  of  +zz  constant  will  demonstrate  a technique 

for  simplifying  the  solution  of  the  three-dimensional  problem.  The 
equation  that  results  from  this  linearization  is 


X+  ♦ ♦ ♦ 

Tzz  Trr 


♦r/r  4 Wr2  " 0 * 


J.  F.  R.  Bailey,  "Nwterioal  Calculation  of  Transom o Floio  About  Slender 
Bodies  of  Revolution, " NASA  TN-D-6S82,  December  1971. 
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where  X is  the  constant  coefficient.  This  equation  is  of  the  form  that 
would  be  obtained  from  a strictly  subsonic  or  supersonic  problem  in 
which  the  Mach  number  is  well  away  from  one.  It  is  not  necessary  to 
resort  to  numerical  techniques  to  solve  this  linearized  equation. 
Ordinary  analysis  may  be  applied.  As  the  problem  is  periodic  in  6 a 
Fourier  transformation  is  a logical  step.  Such  a transformation  results 
when  the  potential  is  expanded  in  a cosine  series  in  0,  as  given  by, 


♦ (z,r,0)  4,5°  (i.r)  ♦ 0 (z.r)  cos  (0)  ♦ C2  (*,r)  cos  (20)  ♦ ... 


A two-dimensional  differential  equation  may  then  be  obtained  for  each 

coefficient  £ . Fortunately,  because  of  the  boundary  conditions,  only 
two  of  these  coefficients  are  non-zero.  The  net  result  is  the  trans- 
formation of  the  original  linear  three-dimensional  problem  into  a pair 
of  two-dimensional  problems,  given  by, 


X£°  ♦ C°  ♦ £°/r  - 0 , 

zz  rr  r * 


surface 


■ dR/dz  , 


£ - o 


and 

UiT  ♦ ♦ Ojr  - cVr2  - 0 , 

zz  rr  r 

C].|  ■ 0 , ■ a r . 

I surface 

As  the  three  dimensional  problem  has  been  transformed  into  two  two- 
dimensional  problems  it  will  be  about  twice  as  difficult  to  solve  the 
three-dimensional  problem  as  it  is  to  solve  a two-dimensional  problem. 
This  represents  a consi  orable  improvement  over  the  factor  of  30 
mentioned  in  the  Introduction. 

Reyhner1*  has  published  a three-dimensional  solution  of  the  full 
potential  equation.  In  his  paper  he  observes  that  the  solution  is 
closely  approximated  by  a constant  plus  a cosine  in  the  circumferential 
direction.  This  result  is  not  surprising  in  light  of  the  above  dis- 
cussion. The  solution  is  exactly  a constant  plus  a cosine  in  the  cir- 
cumferential direction  for  the  subsonic  and  for  the  supersonic 
problems  and  is  not  far  different  for  the  transonic  non-linear  problem. 


4.  T.  A.  Reyhner,  "Transonic  Potential  Flout  Around  Axisyrmetria  Inlets 
and  Bodies  at  Angle  of  Attaokt"  AIAA  Journal , Vol.  IS,  No.  9, 
September  1977 , pp.  1299-1306. 
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The  ease  with  which  the  linear  problem  may  be  solved  in  three 
dimensions  gives  one  hope  that  the  transonic  problem  may  be  similarly 
transformed.  Fourier  transform  techniques  cannot  be  carried  over 
directly  to  the  non-linear  transonic  equation,  however.  A first  require 
ment  for  the  use  of  Fourier  transform  techniques  is  that  the  equation 
to  which  they  are  applied  must  be  linear.  The  transformation  of  the 
non-linear  term  of  equation  (lj  would  yield  a convolution  that  would 
couple  the  resulting  set  of  differential  equations.  It  is  also 
important  to  note  that  a second  requirement  for  the  use  of  Fourier 
transforms  is  that  the  operator  to  which  they  are  applied  must  be 
independent  of  the  choice  of  where  the  6 = 0 plane  is  located.  That 
is,  the  operator  must  have  complete  rotational  symmetry  in  9. 

Fourier  transform  techniques  can  be  applied  indirectly  to  the 
solution  of  the  non-linear  problem.  An  iterative  approach  may  be 
applied  numerically  to  the  solution  of  equation  (1).  In  this  approach 
the  derivatives  with  respect  to  z which  occur  in  the  non-linear  term 
are  evaluated  from  an  old  iteration  of  the  solution.  The  remaining 
two-dimensional  system  in  r and  0 is  linear  and  is  solved  to  obtain  the 
potential  at  the  next  iteration  of  the  solution  on  a plane  in  r and  0. 
The  solution  may  be  obtained  in  this  manner  for  all  of  the  planes 
perpendicular  to  the  z axis  in  the  flow  field.  The  process  may  then  be 
continued  through  many  such  iterations,  evaluating  the  z derivatives 
from  the  solution  at  the  old  iteration  and  obtaining  the  solution  at 
the  new  iteration  from  the  resulting  linear  equations  in  r and  0.  Such 
a process,  of  course,  may  or  may  not  converge.  The  linear  problems 
generated,  however,  are  symmetrical  in  0 and  Fourier  transform 
techniques  may  be  applied. 

This  iterative  technique  in  which  the  potential  is  obtained  for 
planes  cutting  the  body  axis  from  known  values  of  the  z derivative 
terms  is  in  close  harmony  with  the  physical  nature  of  the  flow.  An 
examination  of  the  shadowgraph  shown  in  Figure  1 leaves  one  with  a 
clear  image  of  the  radial  nature  of  the  pattern  of  shocks  about  the 
body.  The  effect  of  a change  in  the  body  geometry  is  felt  far  out  in 
the  flow  field  in  a radial  direction  from  the  cause  of  the  disturbance. 
The  effect  in  the  longitudinal  direction  is  not  nearly  as  great. 

One  would  expect  to  find  this  situation  modeled  in  the  differential 

equation  that  is  used  to  predict  transonic  flow.  The  <j>  derivative 

z z 

term  in  equation  (1),  that  couples  the  potential  at  any  point  to  the 
potential  at  longitudinally  neighboring  points,  is  scaled  by  a small 
coefficient,  1 - M2  - M2  (y  1)  4>z-  When  the  local  Mach  number  is 

one  this  coefficient  is  zero.  Since  the  local  Mach  number  is  always 

near  one  in  transonic  flow,  this  coefficient  is  always  much  less  than 

one.  Near  the  body  surface  the  <f>  term  may  be  neglected  entirely  in 

z z 

regions  away  from  corners  and  other  discontinuities.  The  equation 
that  results  when  the  non-linear  term  is  dropped  is  called  the  inner 
equation  and  is  given  by. 
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(2) 


a ♦ 6 /r  ♦ a „/r2  * 0 
*rr  V veo' 

There  is  no  coupling  in  the  longitudinal  direction  included  in  this 
equation. 

The  radial  propogation  of  physical  influence  helps  to  explain  the 
success  of  line  relaxation  schemes  such  as  the  one  Bailey3  applied  to 
the  solution  of  the  two-dimensional  cylindrical  problem.  The  direction 
of  the  implicit  calculation  in  this  solution  is  maintained  in  the  radial 
direction  which  corresponds  to  the  direction  of  the  major  physical 
coupling.  It  is  also  of  advantage  to  make  use  of  this  physical  coupling 
in  the  three-dimensional  solution  of  equation  (1) . This  coupling  may 
be  taken  advantage  of  if  the  line  relaxation  scheme  that  Baiiey  used 
in  two  dimensions  is  developed  into  a plane  relaxation  scheme  in  three 
dimensions. 

The  implicit  system  of  equations  associated  with  line  relaxation 
is  of  tri -diagonal  form  and  may  be  readily  solved.  The  system  of 
equations  associated  with  plane  relaxation  will  be  penta-diagonal  in 
form  because  of  the  added  coupling  between  the  equations  around  the 
body.  In  such  a system  there  will  be  non-zero  elements  on  the  main 
diagonal  and  on  four  side  diagonals.  Such  a system  cannot  be  solved  as 
directly.  A finite  Fourier  transformation,  however,  will  reduce  the 
three-dimensional  penta-diagonal  system  to  a series  of  two-dimensional 
tri -diagonal  problems. 

There  is  an  additional  advantage  to  the  use  of  Fourier  transforma- 
tions. They  make  natural  the  use  of  spectral  methods.  In  a spectral 
method  the  potential  is  expanded  in  a cosine  series, 

♦ = £°  ♦ C1  cos  (0)  + i2  cos  (20)  ♦ ...  ♦ cos  (k0)  + ... 

The  second  partial  derivative  with  respect  to  0 is  then, 

= - 0 cos  (0)  - 4£2  cos  (20)  - ...  - k2  cos  (k0)  - ... 

Since  the  coefficients  of  cos  (k0)  for  k larger  than  1 are  very  small 
in  the  transonic  problem  it  is  possible  to  evaluate  the  0 derivative 
with  very  few  Fourier  components.  The  necessity  of  keeping  only  a few 
Fourier  components  implies  that  only  a few  grid  points  need  be  main- 
tained in  0,  as  will  be  shown  below.  In  practice  no  increase  in 
accuracy,  of  engineering  significance,  may  be  seen  between  the  solu- 
tions obtained  from  grids  of  4 and  8 points  in  the  0 direction.  In  the 
linear  case  of  low  Mach  number  the  coefficients  for  k larger  than  1 are 
zero  and  only  two  grid  points  are  necessary.  It  would  be  necessary  to 
use  at  least  16  grid  points  to  evaluate  the  second  derivative  of  cos(0) 
with  the  usual  central  difference  formula. 
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The  Fourier  transformations  which  are  useful  in  the  numerical 
solution  of  equation  (1)  are  finite.  The  finite  Fourier  transformation 
of  the  potential  may  be  written  as: 


C = [ E $n  exp  (2irikn/N)]  / . 

n=l 

k n 

The  5 's  are  the  transformed  potentials  associated  with  the  d>  * s given 

at  N equally  spaced  grid  points  on  a circle  around  the  body.  There  are 
as  many  E^'s  as  there  are  <J>n's.  The  potential  may  be  obtained  from  the 

k 

5 's  by  applying  the  inverse  transformation: 


$n  = [ Z £ exp  (-2irikn/N)]  / ^ . 

k=l 


Since  it  is  possible  to  go  backward  and  forward  between  the  trans- 
formed and  untransformed  forms  of  the  potential,  the  finite  Fourier 
transformation  can  neither  add  to  nor  subtract  from  the  information  in 
the  equations.  The  transformation  would  be  described  in  linear  algebra 
as  a change  of  basis.  It  may  be  written  as  a matrix  operator  F whose 
elements  are: 


f,  = exp  (2irikn) 

K|I1 


The  finite  Fourier  transformation,  written  in  matrix  form,  is  then 

K = F * . 

The  transformed  form  of  an  operator  0 such  as  the  one  that  will  be 
obtained  from  the  transonic  equation  is  given  by, 

FOF'1  . 

Fourier  transform  techniques  are  natural  for  the  cylindrical 
problem.  As  only  a few  Fourier  terms  are  needed,  the  transformations 
use  little  computer  time  and  a fast  solution  algorithm  is  possible. 

The  problem  that  creates  significant  difficulty  in  the  use  of  the 
transform  method  arises  in  the  need  to  stabilize  the  iterative  process 
used  to  solve  the  transonic  problem.  Stable  schemes  can  be  developed 
as  will  be  shown. 


IV.  DIFFERENCE  FORMULATION 


The  difference  formulation  used  to  solve  the  three-dimensional 
cylindrical  transonic  problem  utilizes  Fourier  transformations.  The 
algorithm  is  based  on  plane  relaxation  and  a spectral  method  is  used  in 
the  0 direction.  Type  dependent  differencing,  as  developed  by  Murman 
and  Colo5,  is  also  used.  Central  differences  are  used  to  evaluate  the 
4 r derivatives  throughout. 


In  the  following,  A will  be  used  to  imply  a difference  operator  of 
the  type  described  above.  That  is,  Arr$  will  be  used  to  represent  the 

central  difference  form  for  the  second  derivative  of  $ in  the  r 
direction.  A $ will  imply  the  use  of  a central  difference  formula  in 

subsonic  regions  and  the  use  of  a backward  difference  formula  in  super- 
sonic regions.  Further,  Ano$  will  imply  a spectral  evaluation  of  $>n„. 

WU  U tr 

Difference  equations  for  the  transonic  problem  then  may  be  written 
as : 


Arr*  ♦ Ar*/r  ♦ A00*/r2  - -XAz:*  (5) 

where  X represents  the  non-linear  coefficient  1 - M2  - M2  (y  ♦ 1)  A <J>. 

Z 

One  such  equation  will  exist  for  oach  grid  point  in  the  flow  field. 

The  operator  represented  by  tho  left  hand  side  of  this  equation 
meets  all  of  the  criterion  necessary  for  the  successful  application  of 
Fourier  transform  techniques.  It  is  both  linear  and  symmetric  in  e. 

An  iterative  plane  relaxation  procedure  could  be  applied  to  the  solution 
of  the  system  of  equations  of  the  form  shown  in  equation  (5).  This 
procedure  could  be  carried  forward  by  evaluating  the  right  hand  sides 
of  these  equations  from  the  values  of  the  potential  available  from  the 
last  iteration  of  the  solution,  and  by  then  solving  the  resulting 
linear  system  for  the  potential  at  the  next  iteration  of  the  solution. 
Plane  after  plane  of  potentials  could  be  obtained  in  this  manner 
marching  from  a point  upstream  of  the  shell  to  a point  downstream  of 
the  shell  until  the  potential  for  the  entire  flow  field  at  the  new 
iteration  was  known.  Unfortunately,  the  repetition  of  this  process  for 
many  iterations  does  not  converge  to  a stable  solution.  In  order  to 
stabilize  the  iterative  process,  however,  any  term  which  helps  stability 
may  be  added  to  either  side  of  equation  (1)  if  the  term  will  disappear 
when  convergence  is  reached. 


5.  E.  M.  Wumm  and  J.  D.  Cole,  " Calculation  of  Plane  Steady/  Tranaonio 
Floue,"  AIAA  Journal , Vol.  9,  No.  1 , Janutuy  1971 1 pp.  114-131. 
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In  « study  of  convergence  it  is  convenient  to  consider  the  itera- 
tive process  to  be  a marching  process  through  time.  Time  then  forms  a 
fourth  dimension  in  the  problem  and  the  marching  process  may  be  analyzed 
accordingly.  Differences  between  potentials  at  two  consecutive  itera- 
tions are  i hen  considered  to  be  related  to  a derivative  in  time.  On 
convergence,  by  definition,  the  solution  has  reached  a steady  state. 

In  a steady  state  all  time  derivatives  such  as  et  or  mixed  spatial  and 

time  derivatives  such  as  8 will  disappear.  Any  combination  of  9^  or 

0 terms  may  be  added  to  either  side  of  equation  (1)  if  they  will  help 

4 v 

stabilize  the  iterative  procedure.  Ballhaus6  has  discussed  the  effect 
of  0^  terms  as  they  appear  in  the  time  dependent  transonic  small 

disturbance  equation  written  in  Cartesian  coordinates.  In  the  cylin- 
drical problem  it  has  been  found  necessary  to  apply  a type  dependent 
stabilization  process.  The  addition  of  a term  to  the  right  hand 

side  of  the  equation  at  subsonic  points  and  the  addition  of  a r<j>zt 

term  to  the  right  hand  side  of  the  equation  at  supersonic  points  has 
been  found  to  stabilize  the  iterative  process.  The  coefficient  r(r,z) 
is  a suitably  chosen  constant  in  0. 

A stable  scheme  to  which  Fourier  transform  techniques  can  be 
applied  may  then  be  based  on  the  difference  equations  given  by: 

Arr4>  + Ar$/r  + AQ0<ji/r2  = -*Azz4>  + rAxt4>  (supersonic) 

and  (4) 

+ V^r  + A6e^r2  “ "^zz*  + rA-t^  (subsonic) 


In  order  to  see  explicitly  how  convergence  is  achieved  it  is  necessary 
to  write  out  the  right  hand  side  of  equations  (4)  in  detail.  The  nota- 
tion that  will  be  used  is  as  follows.  The  potential  at  the  grid  point 
for  which  the  equation  is  derived  will  be  written  as  <J>.  The  potential 
at  neighboring  points  in  the  z direction  will  be  written  as  or  $ 

for  neighbors  in  the  positive  or  negative  z directions  respectively. 

<p_  will  be  used  to  represent  the  potential  at  the  next  nearest  neighbor 

grid  point  in  the  negative  z direction.  Potentials  which  are  evaluated 
at  the  last  time  step  (iteration)  will  be  marked  with  a A . Equations 
(4)  then  become: 


6.  William  F.  Ballhaus  and.  Howard  Lomax,  " The  Numerical  Simulation 
of  Lou  Frequency  Unsteady  Transonic  Flou  Fields ,"  Lecture  Notes  in 
Physics,  Vol . 35,  pp.  57-63  (1975). 
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Arr*  * Ar*/r  * Aee^r2  * -xU-2$  ♦$_)/a*  ♦ r [ (♦-♦)  - u_  - ♦_)] 

(supersonic) 

and  (S) 

* A A 

Arr*  * Ar*/r  4 Aee^/r2  " 'X^*  ' 2*  * ♦-)/Ar  4 r(*  ‘ *}  * 

(subsonic) 

r(i,r)  - 2Max0(X)/A*  . 


The  non-linear  coefficient  X has  not  been  written  out.  It  is  always 
central  differenced  using  the  old  time  step.  The  coefficient  r is 
dependent  on  the  maximum  value  of  X around  the  circle  of  grid  points  at 
the  current  value  of  r and  2.  Az  gives  grid  separation. 


Equations  (S)  have  been  written  as  they  would  appear  for  any  evenly 
spaced  grid.  The  grid  actually  used  was  not  evenly  spaced.  A slight 
modification  is  necessary  to  write  these  equations  for  a varying  grid. 
This  modification  will  not  be  presented  here.  In  order  to  solve  these 
equations  the  r$  term  is  carried  from  the  right  hand  side  to  the  left 
where  it  becomes  part  of  the  implicit  system  of  equations  which  must 
be  solved.  The  remaining  terms  on  the  right  hand  side  of  equations  (5) 
may  be  evaluated  from  the  solution  for  the  old  iteration  or  from  the 
recently  acquired  solution  on  the  last  plane  for  the  new  iteration. 

The  right  hand  side  is  thus  a known  driving  term  for  the  linear  system 
on  the  left  hand  side.  An  iterative  process,  as  described  above,  may 
then  be  employed. 

Convergence  has  been  found  to  be  monotonic  and  convergence  speed 
may  be  increased  by  over-relaxation,  as  is  used  in  the  line  over- 
relaxation technique.  Speed  of  convergence  may  be  further  increased 
by  use  of  a coarse  longitudinal  grid  for  the  initial  iterations 
followed  by  the  use  of  a fine  grid  when  the  solution  has  partially 
formed.  By  utilixing  these  techniques  quite  reasonable  run  times,  on 
the  order  of  three  minutes,  have  been  achieved  for  three-dimensional 
cases  on  fine  grids. 


V.  DISCUSSION  OF  RESULTS 

The  results  for  computations  of  the  surface  pressure  coefficient 
for  bodies  with  circular  arc  profiles  can  be  seen  in  Figures  3 and  4. 
Figure  3 shows  a comparison  of  computed  and  wind  tunnel  pressure  coef- 


7.  R.  A.  Taylor  and  J.  B.  McDevitt,  "Pressure  Distribution  at  Transonic 
Speeds  for  Parabolic-Arc  Bodies  of  Revolution  Having  Fineness  Ratios 
of  10,  12,  and  14,"  NACA  TN-4234,  March  1958. 


. ■ 


— 


16 


ficient  along  the  surface  of  a 1/10  fineness  ratio  body  at  zero  angle 
of  attack  in  a Mach  number  .99  free  stream.  The  location  of  a shock 
can  be  clearly  seen  just  aft  of  the  center  of  the  body. 

The  solid  line  shows  the  results  of  computations  for  a body 
generated  by  a perfect  circular  arc.  The  wind  tunnel  model,  however, 
was  supported  from  the  rear  by  a sting.  The  effect  of  the  sting  was 
modeled  by  attaching  a solid  cylinder  to  the  body  at  the  location  of 
the  sting.  The  comer  at  the  juncture  of  the  body  and  this  cylinder 
was  faired  to  reduce  the  large  pressure  spike  that  would  have  formed 
at  a sharp  comer.  The  resulting  computed  pressure  coefficients  are 
shown  by  the  dashed  line  in  this  figure.  As  the  angle  of  attack  is 
zero  in  the  case  shown  in  Figure  3,  the  computation  is  two-dimensional. 
This  same  case  was  computed  by  Bailey  in  his  earlier  two-dimensional 
work3  and  the  results  are  identical. 

Figure  4 shows  a comparison  of  computed  and  wind  tunnel8  pressures 
for  a slightly  more  slender  body  of  fineness  ratio  1/12.  The  Mach 
number  in  this  case  was  .90  which  is  too  low  to  allow  development  of  a 
large  supersonic  region  with  strong  shocks.  The  figure  is  presented  to 
show  the  result  of  a three-dimensional  computation.  For  the  results 
shown,  the  sting  was  modeled  in  the  same  manner  as  discussed  above. 

The  results  presented  in  these  two  figures  confirm  the  ability 
of  the  three-dimensional  code  to  predict  surface  pressures  over  smooth 
bodies.  There  is  little  difference  between  the  nose  of  a typical 
artillery  shell  which  is  an  ogive  and  the  front  portion  of  these 
circular  arc  bodies.  Artillery  shell,  however,  often  exhibit  comers, 
particularly  at  the  junction  between  the  ogive  and  cylinder  portions 
and  between  the  cylinder  portion  and  the  boattail.  Strong  shocks  are 
formed  by  the  collapse  of  supersonic  regions  which  are  generated  by  the 
expansion  of  the  flow  over  these  comers  when  the  shell  is  flown  at  a 
slightly  subsonic  velocity  (.8  < M < 1).  As  discussed  in  the  intro- 
duction, the  flow  pattern  generated  by  the  comer  at  the  beginning  of 
the  boattail  is  largely  responsible  for  the  critical  behavior  of  the 
overturning  moment.  Thus,  the  accurate  treatment  of  comer  flow  is  of 
prime  consideration. 

The  ability  of  the  present  theory  to  predict  flow  over  a comer 
can  be  seen  in  Figure  5.  Figure  5 shows  a comparison  of  computed 
surface  pressures  to  wind  tunnel9  experiments  for  flow  over  a cone 
cylinder  model  at  Mach  number  1.1.  The  theory  shows  reasonable 


8.  J.  B.  MoDevitt  and  R.  A.  Taylor , "Force  and  Pressure  Measurements 
at  Transonic  Speeds  for  Several  Bodies  Having  Elliptical  Cross 
Sections , " MCA  TN-4362 , September  1958. 

9.  W.  A.  Page3  "Experimental  Study  of  the  Equivalence  of  Transonic 
Flow  About  Slender  Cone-Cylinders  of  Circular  and  Elliptic  Cross 
Section MCA  TN-4233 , April  1958. 
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behavior  near  the  corner  of  the  cone  and  cylinder  sections.  In  order 
to  achieve  these  results  it  was  necessary  to  use  care  in  applying 
boundary  conditions  at  the  body  surface.  An  approach  that  is  often 
taken  in  the  application  of  boundary  conditions  is  to  use  solutions  of 
the  simpler  inner  equation  (2)  to  extrapolate  the  boundary  conditions 
from  the  body  surface  to  the  body  axis  or  to  some  other  convenient 
location.  In  Bailey's  two  dimensional  paper3  the  boundary  conditions 
were  extrapolated  to  the  axis  where  a series  of  sources  and  sinks  were 
placed.  The  source  and  sink  strengths  were  obtained  from  the  solution 
to  the  inner  equation  (2). 

This  procedure  is  not  feasible  if  accurate  corner  flow  is  to  be 
obtained.  The  equation,  which  is  obtained  by  dropping  the  non-linear 
term  from  equation  (1),  does  not  apply  near  corners  where,  in  fact,  the 
non-linear  term  may  be  large  close  to  the  body  surface.  Boundary 
conditions  must  be  applied  directly  at  the  surface  without  extrapola- 
tion. 


Additional  improvement  in  the  application  of  boundary  conditions 
for  flow  over  corners  may  also  be  obtained  as  indicated  below.  The 
usual  boundary  condition  which  is  applied  at  the  body  surface  is  given 
by, 


♦ 


r 


* dR/dz 

surface 


where  the  left  hand  side  is  the  radial  derivative  of  the  potential 
evaluated  at  the  body  surface  and  the  right  hand  side  is  the  slope  of 
the  body  surface.  This  is  a first  order  approximation  to  the  body 
surface  boundary  condition.  A second  order  formula  is  more  appropriate 
and  is  given  by, 


*r 


surface 


(1  - ♦ , 


) dR/dz  . 

surface 


The  first  order  formula  works  well  as  long  as  2 remains  small  in 

comparison  to  1.  Near  a corner  may  become  large  enough  that  it 

produces  a noticeable  effect  as  seen  in  Figure  S.  Because  of  the 
iterative  relaxation  procedure  used  in  solving  the  potential  equation 
may  be  obtained  at  an  old  iteration.  The  right  hand  side  of  the 

second  order  formula  may  thus  be  evaluated.  The  effect  of  using  a 
finer  grid  spacing  may  also  be  seen  in  Figure  5.  The  grid  in  all  cases 
shown  in  this  figure  was  made  up  of  64  points.  In  the  case  labeled 
fine  grid,  these  points  were  clustered  so  as  to  give  twice  the  density 
near  the  corner.  Subsequent  calculations  have  been  carried  out  with 
128  points  so  as  to  achieve  this  same  fine  density  when  more  than  one 
corner  is  present  on  the  body. 
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It  has  been  shown  in  the  above  discussion  that  the  transonic 
techniques  that  have  been  developed  will  predict  flow  over  both  smooth 
bodies  and  bodies  with  corners.  It  should,  therefore,  be  possible  to 
obtain  the  solution  over  a body  that  closely  resembles  an  artillery 
shell.  The  accuracy  with  which  the  shock  locations  may  be  predicted 
for  a projectile  can  be  seen  in  Figure  6.  The  plot  shown  in  this 
figure  gives  the  pressure  distribution  along  the  surface  of  a projectile 
like  body.  This  body  differs  from  the  M549  projectile  shape  in  that 
it  has  a sharp  nose  and  no  rotating  band.  The  outline  of  an  actual 
M549  projectile  appears  at  the  bottom  of  Figure  6.  The  expansion 
locations  and  shock  locations  are  marked  with  dashed  and  solid  lines 
respectively  about  this  outline  as  taken  from  a shadowgraph10  of  the 
projectile  in  flight  at  M ■ .91.  There  is  a double  shock  pattern 
associated  with  the  collapse  of  each  supersonic  region  as  may  be  seen 
along  the  outline.  The  second  shock  is  due  to  the  interaction  of  the 
first  shock  and  the  viscous  boundary  layer.  The  boundary  layer  separa- 
tes at  the  base  of  the  first  shock  and  reattaches  at  the  base  of  the 
second.  Reasonable  agreement  between  the  shock  locations  seen  in  the 
computed  pressure  distribution  and  the  primary  shock  locations  taken 
from  the  shadowgraph  are  seen. 

The  secondary  shock  is  an  effect  caused  by  the  viscous  boundary 
layer  which  is  not  modeled  in  the  inviscid  computations.  The  presence 
of  the  rotating  band  and  the  effect  of  the  viscous  boundary  layer  in 
rounding  the  corner  at  the  boattail  are  also  not  modeled  in  the  compu- 
tations. Aerodynamic  coefficients  computed  for  the  projectile  like 
body  are  not,  therefore,  oxpected  to  be  in  exact  agreement  with  the 
aerodynamic  coefficients  for  the  actual  projectile.  The  transonic 
effects  can,  however,  be  computed  and  qualitatively  correct  behavior 
can  be  predicted.  The  lift  loading  computed  for  the  projectile  shape 
is  plotted  in  Figure  7.  This  graph  shows  the  normal  force  per  unit 
length  plotted  as  a function  of  the  position  along  the  shell.  It  is 
felt  that  the  features  of  this  curve,  particularly  the  large  downward 
spike  in  the  boattail  region,  give  an  accurate  representation  of  the 
aerodynamic  forces  on  this  body.  Comparison  of  the  pitching  moment 
coefficient  computed  for  this  body  and  range  measurements11  of  the 
pitching  moment  for  the  M549  projectile  are  given  as  a function  of 
Mach  number  in  Figure  8.  The  peak  shown  in  the  computed  results  falls 
a few  hundredths  of  a Mach  number  higher  than  the  peak  in  the  range 
measurements  and  is  not  as  pronounced. 

The  trend  of  design  alterations  may  be  predicted  by  the  present 
technique.  The  boattail  on  the  projectile  shape  discussed  above  was 


10.  Leonard  C.  MaaAllister,  U.S.  Amy  Ballistics  Reseai'oh  Laboratory, 
Aberdeen  Proving  Ground,  Maryland,  shadocAiraph  taken  in  BRL 
Trane onic  Range. 

11.  Robert  L.  McsCoy,  U.S.  Amj  Ballistics  Research  Lal'oratory,  Aberdeen 
Proving  Ground,  Maryland,  private  csornmnioation. 
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shortened  producing  the  results  seen  in  Figure  9.  The  flow  over  the 
boattail  becomes  fully  supersonic  at  a lower  Mach  number  for  the 
shortened  tail  and  transonic  effects  due  to  the  shocks  on  the  boattail 
thus  disappear  at  a lower  Mach  number.  This  causes  the  peak  in  the 
pitching  moment  to  appear  at  a lower  Mach  number. 


VI . SUWARY 

In  order  to  compute  the  aerodynamic  coefficients  for  artillery 
projectiles  it  has  been  necessary  to  resolve  the  shock  patterns 
associated  with  comer  flows.  This  resolution  is  obtained  by  using 
very  fine  grids  and  by  using  care  in  the  application  of  boundary 
conditions.  A practical  computational  algorithm  has  been  developed 
which  is  based  on  Fourier  transformation  techniques  and  which  has  been 
used  to  solve  the  transonic  small  disturbance  equation  on  a cylindrical 
coordinate  system  in  three  dimensions.  The  ability  of  these  methods 
to  accurately  predict  transonic  flow  about  clean  aerodynamic  shapes 
has  been  demonstrated  by  comparisons  of  the  computations  to  experi- 
mental data  for  circular  arc  and  cone  cylinder  bodies.  The  techniques 
have  also  been  shown  to  predict  the  correct  qualitative  trends  of  Cma 

vs  Mach  number  for  actual  shell  configurations.  In  order  to  achieve 
improved  agreement  between  computed  and  experimental  aerodynamic 
coefficients  for  actual  shell,  it  will  be  necessary  to  include  effects 
of  the  viscous  boundary  layer  and  to  include  a model  for  the  rotating 
band. 
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Figure  3.  Comparison  of  Calculated  Pressure  Coefficients  With  Wind 
Tunnel  Data  for  a Fineness  Ratio  1/10  Circular  Arc  Body, 
M • 0.99 
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Figure  4.  Comparison  of  Calculated  Pressure  Coefficients  With  Wind 
Tunnel  Data  for  a Fineness  Ratio  1/12  Circular  Arc  Body 
at  Angle  of  Attack,  a » 4°,  M - 0.90 
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Figure  5.  Comparison  of  Calculated  Pressure  Coefficient  With  Wind 
Tunnel  Data  for  a 7°  Half  Angle  Cone  Cylinder,  M = 0.99 
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Figure  6.  Calculated  Surface  Pressure  Coefficient 
on  an  M549  Projectile  Shape 
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LIST  OF  SYMBOLS 


Cjna  pitching  moment  coefficient  about  center  of  mass 

0^  pressure  coefficient 

f,  matrix  element  of  F 

K,n 

F Fourier  transformation  matrix 

k index  associated  with  Fourier  components  (wave  number) 

M Mach  number 

n index  associated  with  grid  location  in  0 direction 

N total  number  of  grid  points  in  0 direction 

0 arbitrary  matrix  operator 

r radial  coordinate 

R radius  of  body 

z longitudinal  coordinate 

a angle  of  attack 

Y ratio  of  specific  heats  (1.4) 

r coefficient  of  time  derivative  term 

A difference  operator 

AN  lift  loading 

0 circumferential  angular  coordinate 

A non-linear  coefficient 

v transformed  potential 

t 

£ k'th  Fourier  component 

<p  perturbation  velocity  potential 
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